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The FLUKA Monte Carlo program is used to calculate the effects of hadroproduction by primary 
gamma rays incident upon the earth's atmosphere; for the results presented in this paper, only 
primary angles at degrees from zenith are considered. The FLUKA code is believed to be quite 
accurate in reproducing experimental photon hadroproduction data in the 1 GeV to 10 TeV energy 
range studied. The charged pions which are so produced can decay to muons with sufficient energy 
to reach ground level. The number of these muons and their radial and energy distribution are 
studied for incident gamma ray energies from 1 GeV to 10 TeV. The number of these muons is not 
negligible; they can, in certain circumstances, be used to study potential sources of gamma rays like 
gamma ray bursts. It is found, for example, that a 10 TeV incident primary gamma ray produces, on 
average, 3.4 muons which reach ground level; the gamma ray energy which produces the maximum 
number of muons at ground level depends on the spectral index of the primary gamma spectrum, 
a constant which describes how the primary gamma flux rises with decreasing primary energy. An 
example: for a differential spectral index of 2.7, there is a broad maximum number of muons coming 
from » 30 GeV primary gamma ray energy. 

PAC numbers: 95.75.Pq,98.70.Sa,98.70.Rz,13.60.Le 



I. INTRODUCTION 



Primary cosmic rays consist mainly of an isotropic flux of charged particles (primarily protons and nuclei) . Before 
reaching the earth, galactic magnetic fields deflect their paths so information about the angular position of their 
source is lost. On the contrary, neutral particles such as gammas and neutrinos (neutrons decay before reaching the 
earth) can be directly used to locate the angular position of their origin. Therefore gamma astronomy is important 
in the study of well localized exotic astrophysical objects. For example, recent reviews are 0,^]. 

Cosmic gamma rays can be detected directly only by satellite or balloon experiments located essentially outside 
the earth's atmosphere. The detector size and weight sets limits to the sensitivity and energy range which these 
experiments can cover which, at present, extends to several tens of GeV. Higher energies are better investigated by 
means of ground-based experiments which sample the numerous secondary particles produced by the high energy 
primary photons when they interact in the atmosphere. 

Photons constitute only a small fraction of primary cosmic rays. Hadronic showers at ground level are very similar 
to electromagnetic showers because, at each generation of hadronic pion production, about one-third of the pions 
are ir° which immediately decay to gamma rays which then initiate electromagnetic sub-showers. By the time the 
hadronically initiated shower reaches ground level it is essentially all electromagnetic because of the many generations, 
each feeding one-third of their energy into the electromagnetic sector. Therefore, gamma showers at ground level are 
qualitatively similar to those produced by protons and nuclei ||, and experiments based on earth which search for 
gamma primaries must subtract a large background due to the more numerous cosmic protons and nuclei. Various 
methods have been developed to discriminate between the two types of showers |Q,||^J|] . One of the most commonly 
used techniques is based on their different muon content Q . 

Muons are created in the atmosphere mainly as decay products of charged mesons, which are abundant in hadronic 
cascades. Pions and kaons can also be produced in the nuclear interactions of high energy primary gamma rays which 
then initiate secondary hadron showers. The photon cross section for hadronic interactions is about two orders of 
magnitude smaller than that for producing electron pairs. This low hadronic cross section for gamma rays has been 
used as a signature for showers initiated by hadrons. 

However, muons are present in gamma showers too, albeit in smaller numbers than in a hadronically initiated 
shower of the same primary energy. The first estimates confirmed the "muon-poor" characteristics of gamma showers 
However, in 1983, Samorski and Stamm Q reported the results of an extensive air shower experiment at Kiel 
claiming an excess of events with energies above 2000 TeV centered at the angular location of Cygnus X-3. They also 
reported a non-deficiency of muons in these data, seemingly inconsistent with a primary gamma ray hypothesis. In 
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order to evaluate and discuss these data, many authors performed analytical or Monte Carlo calculations of the muon 
flux produced in gamma showers |10|-|2lf . Most of these calculations were one-dimensional and all of them referred to 
gamma energies much larger than 10 TeV, appropriate to the Kiel experiment. The energy realm of these calculations 
was considerably larger than the 1 GeV to 10 TeV region considered in this paper. Therefore, their stress was mostly 
on sources of muons at very high energies (Bethe-Heitler /i-pair production, decay of charmed particles), with less 
emphases on muon decay, energy loss, and radial distribution which become important in the energy range considered 
in this paper. 

More recently, the same calculation techniques have been refined and used to optimize current experiments. In 
addition to an improvement of the adopted cross sections , muon energy loss and decay have been added |^3],^4| , 
in order to better estimate the sensitivity of modern experiments having a lower muon energy threshold. Several 



experiments, GRAND |2{| and Milagro |26| for example, have the capability of measuring the angles of muons at 
ground level with high statistics. With the advantage of large numbers of muons, it becomes possible for these 
experiments to study point sources of gamma rays by studying accumulations of muons at specific angular positions 
p7| |3Cfl . The ability of these experiments to detect a localized source of gamma rays by detecting the secondary muons 
depends upon several factors such as a) the strength, location, duration, and energy spectrum of the source; b) the 
angular resolution, detection area and duration of the experiment; and c) the number of muons which reach detection 
level for the relevant region of primary gamma ray energies due to the physics of the air blanket covering the surface 
of the earth. A precise calculation of (c) in the 1 GeV to 10 TeV primary gamma ray energy region is the subject of 
this paper. 

A reliable study of muon production by primary gamma rays is necessary in these cases to gauge the sensitivity of 
the experiments to gamma primaries and to provide information about the expected energy and spatial distribution 
of muons. Since the threshold energy for detecting the muons has been reduced in these experiments, lower muon 
energies and lower primary gamma ray energies than previously investigated are studied in this paper. Although the 
total number of muons per gamma ray at ground level decreases with lower gamma energies, electromagnetic showers 
produced by photons with lower energies can produce more ground-level muons in total due to the steep increase of 
gamma flux at lower energies (a spectrum dN/dE = E -7 , with 7 = 2.4 ± 0.4 0]) (assuming the primary gamma 
energy is well above photopion production threshold). Therefore, the new calculations presented here include the 
primary gamma ray energy range of 1 GeV < E 7 < 10 TeV and secondary muon or electron energies above 3 MeV. 
For the results presented in this paper, only primary angles at degrees from zenith are considered. 

II. THE FLUKA PROGRAM 
A. Physics 



FLUKA, unlike most Monte Carlo codes used in cosmic ray research |3j|-|40J is not specialized for this particular 
field, but is a multipurpose particle transport program with applications as diverse as proton and electron accelerator 
shielding, calorimetry, medical physics, beam design, high and low energy dosimetry, isotope production, etc. p2| , ^3[ . 
Recently however, it also has been used successfully in space and cosmic ray studies |4l]-(43||. 

In FLUKA, different physical models, or event generators, are responsible for the various aspects (particle type, 
multiplicity, energy and angle) of particle production at different energies. These theoretical models have been directly 
tested against a large amount of nuclear experimental data, and have also been indirectly validated b y c omparisons 
with shower measurements, obtained both at accelerators |44]-^7| and in cosmic ray experiments |^l| , ^8| , ^9| . In partic- 
ular, FLUKA has been shown to predict hadron-generated muon spectra at different heights in the atmosphere with 
good accuracy |5Cfl . 

For the present calculations, the following models are relevant (more details can be found in ]45|]): 




• Hadronic interactions above 4 GeV are simulated according to the Dual Parton Mode l ]5lf| . A list of improve- 
ments to the original Monte Carlo version of the model by Ranft |5^] can be found in |32 

• A cascade pre-equilibrium model is used for hadronic interactions below 3 GeV. The model includes pion and 
kaon production. Between 3 and 4 GeV, inelastic hadron collisions are treated according to a resonance-decay 
model. 

• The Vector Meson Dominance model is used for photonuclear interactions at energies larger than 4 GeV. The 
total cross section is based on experimental photon-proton and photon-neutron cross sections up to and including 
HERA energies and scaled to photon- nucleus interaction according to Bauer et al. |33| ■ Shadowing corrections 
are based on experimental data. The interaction of vector mesons with the nucleus is handled by the Dual 
Parton Model. 
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• Delta Resonance excitation (in the framework of the pre-equilibrium model) is used for photonuclear interactions 
below 4 GeV. 

To illustrate the critical role played by event generators in predicting the muon content of showers, the Feynman-x 
distribution of charged pions (in the lab frame) as calculated by FLUKA is reported in Fig. |l| for both gamma and 
protons at 100 GeV. This figure emphasizes a basic difference between gamma and proton-induced showers: the 
gamma primaries have a larger fraction of high-x secondaries than the proton primaries. 
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FIG. 1. The Feynman-x distribution of charged pions (in the lab frame) for 100 GeV gamma ray primaries compared with 
proton primaries of the same energy as contained in the FLUKA Monte Carlo generator for the primary interaction in air. 



If we disregard the lower total cross section for gamma rays to produce mesons (about a factor 100), then gammas 
appear to be more efficient at producing energetic, forward directed pions, and therefore penetrating muons. The 
gamma primaries give rise to higher energy secondary mesons (on average) and will thus retain a greater fraction of 
the primary energy in the numerous subsequent hadronic interactions for those which do not decay. Although these 
differences seem rather substantive, the actual differences which can be observed near sea level may be minor. At 
any rate, in this paper we only investigate the gamma primaries and leave the comparison with hadronic primaries 
for future consideration. 

The simulation of the electromagnetic cascade in FLUKA is very accurate, including the Landau-Pomeranchuk- 
Migdal effect and a special treatment of the "tip" of the bremsstrahlung spectrum. Electron pairs and bremsstrahlung 
are sampled from the proper double differential energy-angular distributions improving the common practice of using 
average angles. In a similar way, the three-dimensional shape of the hadronic cascades is reproduced in detail by a 
rigorous sampling of correlated energy and angles in decay, scattering, and multiple Coulomb scattering. 

Bethe-Heitler muon pair production is presently being implemented in FLUKA Q,^5), but was not yet available at 
the time of the present calculations; however there is general agreement that this effect is important only for gamma 
energies greater than several TeV and for muon energies larger than 1 TcV pO| , pT| , ^6|Jl^ , |2^ ] . Therefore the results 
presented here should not be affected, with the possible exception of the highest energy point. Charm photoproduction, 
another possible source of muons, is also not available in FLUKA. According to Berezinskii |jq] the corresponding 
cross section at any energy does not exceed a few percent of other muon producing effects. Each of these effects, if 
included, would slightly increase the number of muons at ground level above the numbers obtained in this paper. 



3 



B. Variance reduction techniques 



Statistical techniques for accelerating convergence of the results have been used in a few cosmic ray codes, for 
example MOCCA |3{|. In FLUKA, as in most Monte Carlo codes used to solve deep penetration shielding problems, 
several so-called "biasing" options are available which allow sampling of events having a very small probability. 
Rigorous proofs of the convergence of results obtained by these techniques to the correct value can be found in 
specialized books |]57| , ^8| . It is important to remember, however, that their use is restricted to the estimation of 
expectation values and is not appropriate when studying correlations and fluctuations. 

In this study, the use of variance reduction techniques has proved essential. It is important to realize that the goal 
was not to obtain the same results using less computer time, but to include in the study phase space regions which 
would otherwise not be accessible to Monte Carlo techniques. Due to the very large number of primary photons, some 
interactions, although extremely rare, may generate events having a finite probability to be detected in an experiment. 
Below some level of probability per primary photon the computing time required to collect a sufficient number of such 
events by an unbiased simulation would become prohibitive. 

In the present calculations, the following biasing options have been used: 

Leading particle biasing: At each electromagnetic interaction with two particles in the final state (bremsstrahlung, 
pair production, etc.) only one of the two particles is followed, with a probability proportional to its energy. Its 
statistical weight is modified so as to conserve total weight. This technique, first introduced by A. Van Ginneken 
[ |59t , is very similar - but not identical - to the so-called "thinning algorithm" of Hillas |m| . 

Biasing of the mean free path: For interactions having a very small cross section (in our case photonuclear 
interactions) the cross section is artificially increased by an arbitrary factor chosen by the user (ranging from 
10 to 50 in this calculation). The weight of the produced secondaries is reduced so as to conserve probabilities. 

Forced decay: A similar technique is used to enhance muon production by artificially decreasing the average decay 
length of charged mesons. Also in this case, since the weights of both the parent meson and of the produced 
muon are adjusted by the ratio between the actual and the artificial probability, all the resulting space, energy 
and angular distributions are correctly reproduced (but with much better statistics in the ranges of interest). 

Importance splitting: The loss of statistics due to decrease of particle number with depth in the atmosphere is 
compensated by replacing a Monte Carlo particle with additional identical particles of lower statistical weight 
when the particle crosses a boundary between two regions of different pre-defined statistical importance. 

III. CALCULATIONS 

The 1999 version of FLUKA was used to calculate the muon flux and kinetic energy spectrum at various altitudes. 
At sea level and at 222 m (the GRAND altitude), muon tracks were scored in 10 concentric rings of radii ranging from 
10 m to 10 km. Photons of energy between 1 GeV and 10 TeV were assumed to enter the atmosphere vertically at 
80 km and the full generated electromagnetic and hadronic showers were followed down to pion production threshold 
energy. The mean free path for photonuclear interaction was artificially shortened to enhance the sampling frequency 
by a factor of 10 to 50. As stressed above, the mathematical treatment used ensures that all results are unaffected, 
while the variances of the average scored quantities (due to the statistical nature of the Monte Carlo technique) are 
reduced to acceptable levels. 

Some approximations were adopted which are expected to have a negligible effect on the results. The atmosphere's 
geometry was assumed flat and was subdivided into 50 layers for the muon calculation, each of constant density, in 
order to approximate the exponential character of the earth's air density. Doubling the number of layers from the 
25 used in a previous series of calculations did not show any significant difference. The earth's magnetic field was 
ignored. 

For comparison purposes, a similar set of calculations was made to estimate the electron flux at the same positions. 
All conditions were identical to the muon calculation except the shower energy cutoff was lowered from 150 to 3 MeV 
to accomodate the lower energies of interest for the electrons and the atmosphere was subdivided into only 25 layers. 
The word "electron" in this paper is used generically; in all cases it means e + + e _ ; as well, "muon" means /j + + /i~ . 
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IV. RESULTS 



Fig. ^| shows how the average total number of muons and electrons grows with increasing photon energy. A summary 
of calculated data is reported in Table |. The entries are for gamma ray primaries vertically incident at the top of 
the earth's atmosphere. All error bars in the figures and tables refer to the statistics of the Monte Carlo calculation. 




Photon energy (GeV) 

FIG. 2. Number of muons and electrons at 222 m a.s.l. as a function of incident primary gamma ray energy. The lower 
points are for secondary muons within a 50 m radius. 
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7 energy electrons 
(GeV) 222 m a.s.l. sea level 



muons 

222 m a.s.l. sea level 



10 4 


158 ± 7 


121 ± 6 


3.46 ± 0.05 


3.31 ± 0.05 


3000 


30 ± 3 


23 ± 3 


0.840 ± 0.004 


0.807 ± 0.004 


1000 


4.8 ± 0.7 


3.7 ± 0.6 


0.2312 ± 0.0014 


0.2230 ± 0.0014 


300 


0.51 ± 0.06 


0.41 ± 0.05 


(5.55 ±0.08) x 10" 2 


(5.37 ±0.08) x 10" 2 


100 


(8.79 ±0.06) x 10" 


2 (5.84 ± 0.04) x 10 


~ 2 (1.470 ± 0.008) x 10" 2 (1.424 ± 0.008) x 10" 2 


30 


(1.06 ±0.16) x 10~ 


2 (8.8 ± 1.4) x 10" 


3 (3.25 ± 0.07) x 10" 3 


(3.15 ±0.07) x 10" 3 


10 


(1.4 ±0.3) x 10~ 3 


(1.0 ±0.5) x 10" 


3 (5.66 ± 0.11) x 10" 4 


(5.45 ±0.10) x 10" 4 


3 


(8.7 ± 3.0) x 10 _s 


(3.8 ±1.6) x 10" 


5 (5.42 ±0.09) x 10~ 6 


(4.32 ±0.09) x 10~ 6 


1 






(1.15 ±0.19) x 10~ 9 


(6.9 ± 1.4) x 10" 10 


TABLE I. Number of secondary muons and electrons 


per primary gamma ray within a 10 km radius at 222 and m above 


sea level. 










7 energy 






muons 




(GeV) 


5000 m 


4000 m 


3000 m 2000 m 


1000 m 


10 4 


9.59 ± 0.08 


8.13 ±0.07 


6.56 ± 0.07 5.16 ±0.06 


4.09 ±0.05 


3000 


2.420 ± 0.007 


1.966 ± 0.007 


1.547 ± 0.007 1.218 ±0.006 


0.978 ± 0.005 


1000 


0.6661 ± 0.0019 


0.525 ± 0.002 


0.4105 ± 0.0017 0.325 ±0.002 


0.2659 ± 0.0016 


300 


0.1505 ±0.0014 


0.1175 ±0.0012 


(9.29 ± 0.11) x 10~ 2 (7.54 ± 0.10) x 10" 


" 2 (6.28 ± 0.09) x 10" 2 


100 


(3.727 ±0.017) x 10~ 2 


(2.919 ±0.015) x 10~ 2 


(2.349 ±0.011) x 10~ 2 (1.939 ±0.010) x 10~ 2 (1.647 ± 0.008) x 10" 2 


30 


(7.75 ±0.17) x 10~ 3 


(6.16 ±0.13) x 10" 3 


(5.03 ± 0.10) x 10" 3 (4.25 ± 0.08) x 10" 


" 3 (3.65 ± 0.08) x 10" 3 


10 


(1.46 ±0.03) x 10~ 3 


(1.16 ±0.02) x 10" 3 


(9.4 ±0.2) x 10" 4 (7.71 ±0.15) x 10" 


" 4 (6.43 ± 0.13) x 10" 4 


3 


(1.076 ±0.010) x 10~ 4 


(6.86 ±0.08) x 10" 5 


(4.15 ± 0.05) x 10" 5 (2.31 ± 0.04) x 10" 


5 (1.102 ±0.019) x 10~ 5 


1 


(1.26 ±0.15) x 10~ 6 


(2.2 ±0.3) x 10" 7 


(7±3)xl0~ 8 (5.6 ±0.6) x 10" ! 


5 (2.8 ±0.4) x 10" 9 


TABLE II. Number of secondary muons per incident gamma ray within a 10 km radius at different heights above sea level. 
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Fig. |3| shows how the total number of muons depends upon height above sea level for nine gamma ray primary 
energies and elevations from to 20 km. The corresponding numerical data are contained in Table || for heights from 
1 km to 5 km. One could interpolate among the values in Table || to estimate the expected number of muons per 
gamma primary for ground-based detectors at intermediate heights. As expected, the most probable height for each 
primary energy gradually shifts toward lower heights as the primary energy increases. 




5000 10000 15000 20000 

Height (m) 

FIG. 3. Average number of secondary muons per incident gamma ray within a 10 km radius at different heights above sea 
level. 
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The dependence on radial distance is shown in Fig. ^. The total number of muons within a given radius is plotted 
as a function of radius. For photon energies > 300 GeV, the circle containing half of all the muons reaching sea level 
has a radius < 700 m, increasing to > 1600 m for photon energies < 3 GeV. 
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10 1 10 2 10 3 10 4 

Radial distance (m), R 
FIG. 4. Number of muons per incident photon within the stated radial distance. Nine incident primary gamma ray energies 
are shown. 

A representative example of the numerous spectra which have been calculated is shown in Fig. [5] which shows nine 
different primary photon energies from 1 GeV to 10 TeV and gives the muon energy spectra at 222 m above sea level 
for those muons at distances < 10 km (essentially all of them). For a given primary gamma ray energy the points and 
their error bars are highly correlated due to the summations involved in getting the integrals presented. 
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10 -1 1Q 1Q 1 1Q 2 1Q 3 1Q 4 

Muon energy, E M (GeV) 
FIG. 5. Integral spectra of muon kinetic energy at radial distances < 10 km from the shower center, at 222 m a.s.l. 

To estimate the total number of muons reaching detection level, the calculated data should be folded with the 
incident primary photon spectrum. Assuming a differential energy dependence of E -2 41 J3l| , it is found that the largest 
number of muons comes from primary energies of 27 GeV. If the spectrum were softer (spectral index, 7 = 2.77), 
the corresponding primary energy would drop to 12 GeV; if the spectrum were harder (7 = 2.05), the muon number 
would continue to increase as the primary energy increases and would depend on where the spectrum softens or cuts 
off. 



V. CONCLUSIONS 



These reported results are from the Monte Carlo program FLUKA which contains event generators built on the 
accurate experimental data available in the the energy region of interest. The use of biasing techniques enabled 
the calculation of events originated by low energy primaries with muon probabilities as low as ~ 10 -9 per incident 
gamma ray, which are important since the primary cosmic ray flux rises rapidly with decreasing primary energy. This 
FLUKA-based Monte Carlo calculation is believed to be quite accurate in this energy region. As such, the results 
give a good estimate of the number of muons expected from gamma ray primaries in this energy region and help 
overcome the stigma that muons are only anticoincidence signals for gamma ray primaries. Instead, with the high 
statistics of the experimental results that are becoming available in this 1 GeV to 10 TeV primary energy region, a 
statistically significant excess accumulation of muons at a specific angular direction becomes a positive signature for 
gamma ray initiated hadronic showers and is thus a tool to fill the ener gy gap between satellite/balloon experiments 
and Cerenkov arrays. Examples of this use are contained in references |p0|j60| . 



Thanks to S. Roesler, M. Dunford, and T. Bowen for their help with this paper. Part of this work was supported by 
the Department of Energy under contract DE-AC03-76SF00515. Project GRAND is funded through grants from the 
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